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O Magnetic fields in an accretion disc are advected with the 
accretion flow. The net poloidal flux (flux passing through 
the disc) can not be changed by small scale processes in- 



ABSTRACT 

We investigate the stability to nonaxisymmetric perturbations of an accretion disc in 
which a poloidal magnetic field provides part of the radial support against gravity. Inter- 
change instability due to radial gradients in the magnetic field are strongly stabilized by 
the shear fiow in the disc. For smooth field distributions this instability is restricted to 
discs in which the magnetic energy is comparable to the gravitational energy. An incom- 
pressible model for the instability akin to the Boussinesq approximation for convection 
is given which predicts the behaviour of the instability accurately. Global axisymmetric 
disturbances are also considered and found to be stable for a certain class of models. 
The results indicate that accretion discs may be able to support poloidal fields which 
are strong enough to suppress other forms of magnetic instability. These strong and 
stable field distributions are likely to be well suited for the magnetic acceleration of jets 
and winds. 
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INTRODUCTION 



side the disc itself. This flux is 'inherited' by the disc from 
the accreting matter. Compressed by the accretion flow, this 
flux can build up substantial field strengths in the central 
regions. At least in the case of protostellar discs, the enviro- 
ment in which they form are known to be highly magnetic 
(Heiles et al. 1993), and the question how, while forming a 
star, the accretion flow disposes of the 'excess flux' of mag- 
netic field lines threading the disc is a classical problem in 
star formation (see the review by McKee et al. 1993). 

Depending on their strength and radial distribution, 
poloidal fields can be of the right strength and configu- 
ration to produce magnetically accelerated jets and winds 
(Konigl and Ruden 1993, Blandford 1992). Acceleration of 
these flows works best in strong fields, while their degree 
of coUimation is probably dependent on the radial distribu- 
tion of the poloidal flux (cf. Spruit, 1994). The strength of 
the field is limited by the rate at which the field diffuses 
outward, either by some 'turbulent' process (e.g. Van Bal- 
legooijen, 1989), or by microscopic processes (in the case of 
protostellar discs). Of particular interest is outward diffu- 
sion due to instabilities driven by the magnetic field itself. 



which is expected to become important especially at high 
field strength. Such instabilities are expected, in the form 
of interchanges whenever the magnetic forces provide part 
of the support against gravity. In the same way as in solar 
prominences (Anzer 1969, Priest 1982), energy is gained in 
this case by exchange of the magnetized gas providing the 
support against the less strongly magnetized gas being sup- 
ported. In this form of instability, neighboring field lines with 
the plasma attached to them are interchanged, with little or 
no bending of the field lines. For poloidal fields threading a 
disk, this means that the unstable displacements take place 
in the disk plane, with little dependence on the direction per- 
pendicular to the plane. This interchange process is known 
to be quite effective at transporting mass across field lines 
in solar prominences, where it can be studied in detail (cf 
Priest 1982, Schmieder 1989). The process characteristically 
takes the form of small scale (local) interchanges. It has pre- 
viously been studied in the context of rotating, magnetized 
stars by Pitts and Tayler (1985, and references therein). 

In the context of accretion discs, interchange instability 
has been studied by Spruit and Taam (1989, hereafter ST), 
who concluded that uniform rotation has no effect on the 
instability condition or the maximum growth rates (see also 
Lepeltier and Aly, 1995). From this, they concluded that 
transfer of mass across field lines may well be quite effective 
in the magnetospheres of accreting neutron stars. These re- 
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suits are not applicable to fields in accretion discs, since they 
hold for uniform rotation only, and do not take into account 
the strong shear of the flow in a disc. Shear is likely to have 
a strongly stabilizing effect, since it distorts unstable dis- 
placements, on the short orbital time scale, into less rapidly 
growing ones. 

In this paper, we study the effect of shear on interchange 
instability. We first point out the close analogy between this 
process and 2-dimensional convection, by showing that in 
an approximation akin to the incompressible (Boussinesq) 
approximation in convection, the two problems are mathe- 
matically identical. By integrating the full time dependent 
equations for linear disturbances in a local (shearing sheet) 
model we show that the instability behaves, in the short 
wavelength limit, as expected from an incompressible model. 
The effect of shear turns out to be quite strong, such that 
the instability becomes important only when the magnetic 
support of the disc becomes comparable to gravity, or when 
the field strength varies on a short radial length scale. 

Other forms of instability of poloidal fields are known, 
in particular the rotating-shear flow instability studied in 
the context of discs by Balbus and Hawley (1991). This 
instability is suppressed already when the field strength is 
comparable to the gas pressure. An important conlcusion is 
therefore, that poloidal fields larger than this value are sta- 
ble to both the Balbus-Hawley form and to interchanges, up 
to the point where the field starts contributing significantly 
to the support of the disc. The strong poloidal fields that 
are optimal for the production of jets and winds may thus 
be a natural consequence of the concentration of fields by 
the accretion process. 

A preliminary treatment of the present problem is re- 
ported in Lubow and Spruit (1995). 



2 EQUATIONS 

In the thin disc limit, neither the gas pressure nor the mag- 
netic pressure contributes, and only the magnetic curvature 
force appears in the equation of motion. Since to first ap- 
proximation the interchanges do not bend field lines and are 
nearly parallel to the disk plane, we can take the flow field u 
to be independent of height in the disk. As shown in ST, the 
MHD equations can then be reduced to two dimensions by 
integrating vertically across the disk. In the absence of vis- 
cous forces, these thin-disc MHD equations are, for motions 
in the disk plane: 

du _ [B]B, 
At 

0, (2) 



d(E/5,)/dt = 0. 



(4) 



47rE 
(9tE -|- div(Eu 



+ 9, 



(1) 



dtB, + Ab/{B,u) 



0, 



(3) 



where [E\ denotes the jump in the field vector across the disc 
and g is the acceleration of gravity. The disc is assumed to be 
flat and perpendicular to the z-axis and to remain so during 
perturbations. That is, we do not consider 'buckling modes' 
or corrugation waves (for more general equations allowing 
for corrugations see ST). Due to the thin disc assumption, 
all structure inside the disc has disappeared from view, and 
the curvature force, integrated across the disc appears as the 
first term on the RHS in (1). The continuity and induction 
eqs. (2,3) can be combined into: 



To complete the equations, a prescription is needed to 
determine the field configuration above the disc from the 
values of the normal component on the disc. We assume 
the field above the disc to be a potential field: 



0, 



with boundary conditions 

<9^$ = B^ {z = 0), $ ^ [z ^ oo) 



(5) 



(6) 



In general, one would like to solve for a force-free field above 
the disc, but this is a much more demanding problem. In- 
stead, we assume here that the Alfven speed above the disc 
is sufficiently high that efficient reconnection to a nearly po- 
tential field takes place (as is the case in the solar corona, 
for example). 

Conceptually, the motions determined by the equation 
of motion carry the vertical component of the field around 
by the induction equation (3); the potential field problem 
(5) determines from B^ the horizontal components {Br, B^) 
at the disc surface. These components then determine the 
magnetic force in the equation of motion, thus closing the 
set of equations. 

2.1 Analogy with 2-D hydrodynamics 

The form of the equations has a certain similarity to the 
ordinary hydrodynamic equations in two dimensions. This 
is seen by noting that 

(V X [B]). = 2u = (7) 

by the potential field assumption. Hence we may write 

[B]/47r = -Vh/, (8) 

where Vh denotes the gradient taken in the plane of the 
disc. Writing 



eqs. (1) and (4) can be written as 

du „ r 

"^-TT = - Vh/ -I- rag, 
At 

AmjAt = 0. 



(9) 

(10) 
(11) 



Equation (10) is identical to the Euler equation, with m 
playing the role of the gas density p and / that of the gas 
pressure p. The form of the 'pressure' / is very different 
from the ordinary hydrodynamic case, however. Whereas 
in hydrodynamics one can assume, say, a polytropic rela- 
tion between p and p or an ideal gas, / is related to m in 
a much more complicated way, involving eqs. (11), 3 and 
the solution of the potential problem (5). Hence the anal- 
ogy with hydrodynamics is really close only if additional 
assumptions are made. The problem defined by (1-3), (5) in 
general has compressive solutions (i.e. divu ^ 0), mediated 
by the 'pressure' / (ST, Tagger et al. 1990). The interchange 
instabilities studied in this paper, however, are basically in- 
compressive in nature, and it turns out to be illuminating 
to describe them in an incompressible approximation which 
filters out the stable compressive waves. This can be done 
by taking the curl of (10) to eliminate / and adding the 
approximation divu = 0. In this way, the incompressible 
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limit of the magnetic instability problem becomes equiva- 
lent to 2-dimensional (z-independent) Bousinesq convection. 
The incompressible limit is discussed further in section 3.2. 

2.2 Linearized equations in disc geometry 

We write equations (1,3,4) in polar coordinates (r, (f)) and 
linearize about a basic state Ur = 0, u,f, = Q(r)r. In the 
absence of currugations of the disc, and of externally im- 
posed fields, the field is antisymmetric (dipole-like) about 
the midplane, such that Bz is the same on both sides, while 
the values of Br and B,f, are opposite. Thus we may write 



[Br 



25: 



m = 2Bj 



(12) 



where the index denotes evaluation at the upper surface 
of the disc. We denote by a prime ' the (Eulerian) pertur- 
bations while unprimed quantities refer to the equilibrium 
state, which satisfies 



= gm - g + r, 



where 

ffm = 



B+Bz 
27rE 



(13) 



(14) 



is the magnetic acceleration and g = GM/r^ the gravity due 
to the central object. Due to the magnetic support, Q differs 
from its Keplerian value. The shear rate in this equilibrium 
is 



dQ 3 ffm d 2 

b = r — — = U — — In r 

dr 2 2Qdr ^ " 



(15) 



This also differs from the Keplerian expression 5k 
— 3/2Qk. The linearisation of (1,3,4) yields 

(9t + nd^)u'r - 2Q«; = + 



27rE 27r V E 

(d, + Q9^)«; + u'r-^i^) = 

T At 27rE 

{dt + nd4,)Bz = -div(5,«'). 



(16) 

(17) 
(18) 
(19) 



2.3 Shearing sheet model 

An approximate model for a differentially rotating disc is 
the shearing sheet (Goldreich and Lynden-Bell, 1965). In 
this model local cartesian coordinates (x,y) in the radial 
and azimuthal directions are used around a point r = ro , in 
a frame rotating at the rate Qo = ^{to), while the shear is 
approximated as constant, i.e. 



ro + X, 



l/rod4, dy, 



Q — Qo = Sx I To ■ 



(20) 



This model describes a small section of a disc near r = To, 
such that the curvature of the rotating flow can be neglected, 
and the r-dependence of the equilibrium quantities is as- 
sumed to be linear. In this approximation, the equations 
become 

(9, + 5x9,)«;-2Qo< = ^^^ + ff (^)', (21) 



(dt + Sxdy)u'y + (2Qo + S'K = " "° . 

ZTTZjO 



(22) 



{dt + Sxdy}Bi = -div(5,u'), (23) 

(.. + 5...)(f)' + „;(A|)^ = o, (24) 

where the index o denotes evaluation at r = ro . 

An approximation is needed to simplify the poten- 
tial problem (5). As in ST, we use a short- wavelength ap- 
proximation, in which the length scale of the perturbation 
is assumed to be small compared with the length scale 
of the background gradients. Thus, if the wavenumber is 
k = (kj;, ky ), the field perturbation has the form 

B'(x, y, z, t) = b(t) exTp(ikxX + iky-y — k\z\), (25) 

where 

k = {kl + kl)"\ (26) 
k 

hx = -i-rizSgn(z), (27) 
k 

k 

by = -i-^bzsgn{z). (28) 

We transform to a shearing coordinate frame x' = x, 
y' = y — Sxi, z' = z, t' = i. Hence dt -\- Sxdy = d't, dx = 
d'x ~ st'd'y. In these coordinates eqs. (21-24) have solutions 
of the form 

BzU = BzoU{t') exTp(ikxx' + ikyy'). (29) 
The linearized equations then take the form 



dt'Ux 



dt'Uy 



27rEo k 
B kij 



■l] + F + 2^loUy, 



27rEo k 
dt^lS = -kU 



/? - (2Qo + S)Ux 



phere 



27rEo dr \B 



(30) 

(31) 
(32) 
(33) 

(34) 



The wave vector is now time dependent due to the shearing 
coordinate frame: 



(35) 



which implies fc^ = at t = 0. As long as 5 7^ 0, other 
values of kx{0) are equivalent to a differerent choice of the 
origin t = 0. In the uniformly rotating case 5 = 0, fc^ is a 
constant (generally nonzero) parameter. 

The terms involving /? in the equations of motion rep- 
resent the restoring force due compression of the field lines. 
Since this force is mediated by the external potential field, 
through the curvature force rather than by the magnetic 
pressure, these terms depend in a different way on the 
wavenumber than the normal pressure terms. As a conse- 
quence, the phase speed of waves mediated by these curva- 
ture forces is proportional to k^^^ . 

By introducing shearing coordinates, we have removed 
the linear dependence of the shear flow on x, trading it in 
for an explicit time dependence of the equations. Since the 
waves we study are not exponentially growing, their detailed 
time dependence has to be taken into account anyway (cf. 
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also the discussion in Ryu and Goodman, 1992). Note that 
solutions of the form chosen (in shearing coordinates) re- 
quire that the medium can be regarded as unbounded, or 
that periodic boundary conditions can be adopted in the 
direction of the shear (x-direction) 

From now on, we drop the subscripts o on the quanti- 
ties evaluated at To - The solutions of eqs (30-35) depend on 
only three independent parameters. This is seen by writing 
them in nondimensional form. The system has a character- 
istic length scale: 



L 



Bi 



27rEQ2- 

The ratio i/r can be written as 
Bz gra 



LI 



B+ 



(36) 



(37) 



which measures the strength of the magnetic relative to the 
rotational acceleration in the equilibrium state. Writing 



k = K/L, 
U = QLtl, 



t' = r/Q, 
F = Q^LF, 



the equations become, after dropping the tilde's: 
d,U^ = i^i] + F + 2Uy, 

d,Uy = i^li -{2 + s)U^, 



drl3-- 

drF 



-K U, 



phere 
= 



" Q2 dr Bz 



(38) 
(39) 

(40) 

(41) 

(42) 
(43) 

(44) 
(45) 



The parameters of the problem are thus Ky, a and s. The 
value of Kxit) is fixed, apart from an arbitrary zero point 
in time, by the evolution equation (35). In (44) we have 
introduced, for future reference, the magnetic buoyancy fre- 
quency N„:^. It is the frequency of interchange displacements 
in the case of a stable magnetic gradient, analogous to the 
buoyancy frequency in a convectively stable stratification. 

For the shearing sheet model to be applicable, the wave- 
length of the unstable flow must be shorter than the length 
scale on which the unstable gradient is present, at least dur- 
ing the time that growth takes place. Thus we must require 
that 



dr Bz 



(46) 



In terms of the dimensionless wavenumber K, this implies, 
approximately. 



K^\a\. 



(47) 



the wavenumber does not depend on time in this case, eqs. 
(40-43) have solutions of the form exp((7r), with dispersion 
relation 

2 



(J A + a , a ( K, 

, H {— h 1) H — 

K K ' K V K 



0. 



(48) 



From the sign of the last term it follows that the condition 
for instability is a < (except in the axisymmetric case 
Ky = 0, which is discussed in section 3.3). By consistency 
condition (47) we must have K ^ |a|, hence the last term 
in (48) is small, though K itself can still be large or small 
compared with unity. The solutions are thus approximately 



-(A' + 4), 



2 
(^2 



Ky^' 



K 



K J 4-1- K 

In the high wavenumber limit K 1, this simplifies to 



(49) 



2 



2 
<^2 



(50) 



The first solution is a stable compressive wave, in which the 
restoring force is due to the perturbations in the curvature 
force. In terms of the dimensional quantities, the frequency 
of this wave, in the absence of rotation, is 



'27rE' 



(51) 



Since the restoring force is mediated by the potential field 
outside the disc, the wave frequency has the dependence 
k^^^ which is characteristic also of the very analogous case of 
selfgravitating discs. The second solution is the interchange 
mode, unstable when a < 0, with the growth rate maximiz- 
ing for Ky ^ Kx- It is stable for a > with a behavior 
similar to internal gravity waves. For further discussion of 
the physics of these results, see ST. At wavenumbers K ^ 1, 
where the compressive wave frequency drops below the ro- 
tation frequency, rotation provides an additional restoring 
force in the compressive wave, and reduces the growth rate 
of the interchange mode without, however, changing the sta- 
bility condition. 



3.2 The incompressible limit 

Like in the case of ordinary convection, interchange insta- 
bility is described correctly, at large wavenumber, by the 
incompressible (Boussinesq) limit. As in the case of ordi- 
nary convection, this limit is obtained by taking k ^ oo 
while keeping the 'pressure' force finite [the /9-terms in (40- 
43)]. It has the advantage of showing the instability without 
interference from compressive waves. With K = K this 
limit gives 



k-Ur~.l/K (K ^ oo), 



(52) 



Hence Uy -K^/KyUx- We multiply (40) by Ky, (41) 
by Kj; and subtract (thus effectively taking the curl of the 
equation of motion). Using (52) and (43) to eliminate Uy 
and F then yields 



3 SOLUTIONS 

3.1 Uniform rotation 

The uniformly rotating case has been studied before by ST 
(see also the discussion in Lubow and Spruit, 1995). Since 



— — — -|-4rs — h (2s +a)Uj: 

dr^ dr 



0. 



(53) 



Hence the natural time scale in the incompressible case is 
just the shear time scale 

V = Ts. (54) 
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In this time coordinate, and writing 

f = {l + v')U, 

this equation becomes 



(1 + " )xT + ~f '■ 



0. 



(55) 



(56) 



The equation is second order in time instead of fourth be- 
cause by the incompressibility limit the compressive waves 
have been eliminated. The general solution of (56) can be 
written in terms of hypergeometric functions. Here we need 
mainly the asymptotic dependence of the solutions for large 
time, which are powers of v. The fastest growing solution 
has 



[(l-Aa/s^y^ - 3]/2. 



(57) 



For s = the index n diverges, in accordance with the fact 
that in the uniformly rotating case the growth is exponential 
rather than algebraic. The condition that Ux be a growing 
power of V yields, with the definition of a 



-Nl > 2S\ 



(58) 



The azimuthal velocity grows faster by one power of v, since 
by the incompressibility Uy ~ Kx/KyUx- Thus Uy grows 
in the incompressible model whenever < 0, but only 

linearly or slower if (58) is not satisfied. It is a matter of 
debate whether such slow growth in Uy only should be still 
called unstable. For uniform rotation, (58) reduces to the 
condition —N^ > 0, in agreement with the results in ST. 
With (15), (44) condition (58) translates into 



d , E / dQ\^ 



(59) 



If the length scale on which 'E/Bz varies is of the order 
r, this means that the magnetic support of the disc must 
be comparable to the rotational support, before instability 
sets in. This implies that the magnetic energy of the system 
(disc plus external field) must be comparable to the gravita- 
tional binding energy. Only if the field strength varies with 
r on a smaller length scale does interchange instability take 
place at smaller field strengths. For a disc with finite inner 
and outer radii, this condition also says that the intrinsic 
(nonrotating) instability time scale must be shorter than 
the shear time scale across the whole disc. It is probable 
that additional globally unstable modes would exist in such 
discs with boundaries. 

We conclude that the incompressible limit predicts alge- 
braic growth of interchange modes but only if the magnetic 
field is strong enough to start contributing to the radial sup- 
port of the disc against gravity. 

3.3 Numerical results 

When the assumption of incompressibility is dropped, the 
full 4"* order system (40)-(43) must be integrated. Since the 
instability can be transient in the presence of shear, a growth 
rate or even an algebraic growth index is not necessarily a 
good measure of the effects of the instability. The growth of a 
small disturbance can be limited to a finite time interval, but 
the amplification factor over this interval can be quite large. 
This happens, for example, in the nonaxisymmetric form of 
the instability studied by Balbus and Hawley (1991). 



For this reason, we measure the instability by the ampli- 
fication factor of an initial disturbance over a finite interval 
in time. We determine this factor by integrating the linear 
equations (40-43) in time numerically. We start at t = 
coresponding to Kx = with the initial conditions 



Ux 



1, 



13 = 0, 



F = 0. 



(60) 



In this way the perturbation starts as a purely incompressive 
(K ■ U = 0) motion. Other choices yield a larger admixture 
of stable compressive waves, but very similar amplification 
factors, except at low amplification when the instability is 
rather marginal anyway. In Figure 1 the time development 
is shown, and compared with the incompressible limit dis- 
cussed above. The incompressible results were obtained by 
integrating eq. (53) with the initial conditions 



Ux 



1, 



dUx/dT 



0. 



(61) 



These are the equivalent of (60) for the incompressible case. 
It is seen that the incompressible model is quite accurate 
when Ky is large enough that condition (47) for consistency 
of the shearing sheet approximation is satisfied at all t. The 
asymptotic behaviour is accurately reproduced even when 
(47) is not satisfied initially; this is because for large t, the 
value of K always becomes large enough that (47) holds. The 
comparison indicates additional, transient, growth around 
t = for the smaller azimuthal wave numbers, in the com- 
pressible case. This may be a real effect of compressibility 
(some form of 'swing' amplification), but we can not estab- 
lish this from the present model, since it happens only under 
conditions {K{t) ^ |a|) when the consistency condtion (47) 
is not satisfied. 

We can conclude that the incompressible model seems 
to give the correct description of the instability whenever 
the wavenumber Ky is large enough that the shearing sheet 
model applies. 

Fig 2 shows the amplification factors after integra- 
tion until t = 10, as functions of the driving parameter 
a = N^/fl"^, the shear rate s = S/fl, and the dimension- 
less azimuthal wavenumber Ky. 

The amplification factor shows two different forms of 
instability. At large Ky we expect to find interchange insta- 
bility essentially in its incompressible form, as discussed in 
section 3.2. This is shown in Figure 3a, for Ky = 10. The 
amplification factor show a broad hump centred at uniform 
rotation, s = 0. Amplification factors obtained from inte- 
grating equation (53) which is valid for the incompressible 
limit are shown for comparison in Fig. 3b. The agreement is 
good, in particular the incompressible instability condition 
(58) appears to be valid. 

While this part of the instability behaves as expected, 
the amplification factor shows an additional form of instabil- 
ity at low Ky. An example is shown in Fig. 4, for Ky = 0.25. 
The growth appears to maximize at Ky = 0, the axisymmet- 
ric case. This is seen in Fig. 5, which show the amplification 
factor as a function of wavenumber Ky and instability pa- 
rameter a. 

The axisymmetric case is easily treated analytically 
since in this case the radial wave number Kj; is time in- 
dependent. From eqs (40-43) we have, with Ky = 0, the 
dispersion relation 



-K, 



2{2 + s} 



(62) 
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Figure 1. Velocity amplitude Ux in interchange instability at 
N^/Q^ = —20 and -10, for s = —1.5, showing algebraic growth. 
Time in units of . Dashed: growth predicted by the incom- 
pressible model, top: Ky = 20, bottom: Ky = 200, 



1 O „ 




Figure 2. Amplification factors after Qt = 10 as functions of the 
dimensionless azimuthal wave number, for a = N^/Q^ = —4 — 2s. 

In terms of the physical (dimensional) quantities, this pre- 
dicts instability whenever 

-Nl>2^^{^lr') + kx^. (63) 
r ar /ttL 

The growth rate is of the order 57 except close to marginal 
stability, and instability first sets in at the lowest wavenum- 
bers. The instability is therefore an intrinsically global one if 
it exists. This raises the question whether it can be treated 



2.S ^ 
Z.O ■- 




Figure 3. a: Amplification factor A after fit = 30 for Ky = 10 in 
terms of 7 = In A/30, showing the stabilization of the interchange 
instability away from uniform rotation (s = 0). b: same derived 
analytically in the incompressible limit 



m 




Figure 4. gray scale representation of the amplification factor 
(white highest) for a = —3. 

correctly in the short wave approximation, on which (40-43) 
are based. Comparing (62) with (47) we see that the insta- 
bility is excluded by condition (63) whenever the short wave 
approximation applies. 

The last term in (63) is equal to a;^, where uJc is the fre- 
quency of the compressive wave of wave number kx defined 
by (51). To see the physical interpretation of this term, con- 
sider the radial length scale im over which the destabilizing 
field gradient exists (length scale of iVm). The prospective 
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Figure 5. Amplification factors after t = 5 for s = —3/2 



unstable wave must still fit into this length, for instability 
to occur, i.e. kx > 1/im. Condition (63) then says that for 
axisymmetric instability to occur the intrinsic (nonrotating) 
growth must be faster than both the epicyclic frequency and 
the time for a compressive wave to cross the region with the 
destabilizing field gradient. 

To study the long-wave form of instability properly, a 
more global form of stability analysis is needed. This is, in 
general, a cumbersome problem because of the nonlocal na- 
ture of the magnetic forces. A global result with which cer- 
tain cases can be shown to be stable towards axisymmetric 
perturbations is given in the next section. 



4 AXISYMMETRIC STABILITY 

Analysis of the global stability problem is helped greatly 
by the indication, from the local analysis presented in the 
previous sections, that the low-wavenumber form of insta- 
bility is essentially an axisymmetric one. This indicates that 
its mechanism, if it exists, must be simple. Assuming ax- 
isymmetry from the outset, the additional complication of 
transient instability due to shear disappears. We restrict our- 
selves to discs in which the polarity of Br does not depend 
on r. If the polarity of Bz is also constant, this means that 
the sign of the magnetic acceleration g^a is constant, but it 
is not necessary to assume this. 

The equation of motion (1) can be written (in an inertial 
frame) as 



(64) 



du J X B 

d7 " ~ ^ 

where J(r,4>) is the current integrated over the disc thick- 
ness, and g the gravity due to the central object. In the 
equilibrium state, this yields: 

J^B, _ GM 



(65) 



where the current is related to the radial field component 
by: 



5+/27r 



(66) 



and Br is the radial field at the top surface of the disc, as 
in section 2.3. Because B^ is of constant sign by assump- 
tion, and the stability of the field is unaffected by an overall 



change of sign of the field, we can take the signs of B^ and 
to be positive. Equation (65) defines the equilibrium ro- 
tation rate ^{r). Since we are considering an axisymmetric 
equilibrium with axisymmetric perturbations, the vertical 
component B^ of the field can be written in terms of the 
vector potential A as 

5. = i^, (67) 
T or 

while the vector potential can be written in terms of the 
current distribution in the disc as 



A4, = O(J^), 

where the self-adjoint operator O is defined by 



0{U) = ^ Ids 



cos((/i — 4>')J^{r' , 4>') 



2rr'cos((t> - ,/i')]i/2 



(68) 



(69) 



The integral is taken over the surface s of the disc. The 
operator has the property that 0{F) is positive definite if 
F is positive definite, hence A,f, is positive definite. 

Axisymmetric perturbations are not sheared, and their 
linear development can be described in terms of normal 
modes. Linearizing (64), setting d / dt = 7, the radial and 
azimuthal components can be combined to yield 



2g 



(70) 



where ^ is the displacement (u = d^/dt), a prime denotes 
a perturbed quantity and the epicyclic frequency k is given 
by 



^2 ^ 2Qi^(Qr") 



/,From the induction equation dB/dt 
have, using axisymmetry, 

A'^ = -S.rB, = c(j; 



(71) 

V X (u X B) we 

(72) 



Thus can be written in terms of fr, so that (70) can be 
written as 



(73) 



where C is the self-adjoint operator defined by the RHS of 
(70). Multiplying (70) by E{* and integrating over s we have 



dsElfr 



ds 4*0(4) = -W. 



(74) 



For stability, it is sufficient that the energy W be positive 
definite. The main problem is the estimation of the magnetic 
energy term 



dsj':o{A] 



(75) 



which is positive and hence stabilizing. It corresponds to 
the restoring force that develops (via the potential field 
above the disc) when the field lines threading the disc 
are compressed together or expanded. For nonaxisymmet- 
ric perturbations, the energy expenditure in such compres- 
sions/expansions can be minimized by taking rotational dis- 
placements (i.e. interchanges), but this is not possible for 
axisymmetric perturbations. To find a sufficient condition 
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for stability, we need to find a sufficiently large minimum 
value for (75). To do this, consider the variational problem 
of finding the minimum of /t defined by 



/,ds/(r)|A;p /^ds/(r)|A;p 



(76) 



for some positive definite /(r). The minimum of /t is the 
smallest eigenvalue of the equation 



If we choose 
/(r) = J^IA^ 



(77) 



(78) 



(which is positive definite since and A^ are), this equation 
can be written as 



A^O-^A'^,) = nA(t,'0-\A4,). 



(79) 



It is apparent that (79) has the solution A'^ = A^ with 
/t = f and the corresponding = J^. The current density 
for this eigenfunction is everywhere positive as it must be 
for the smallest eigenvalue /t which is therefore unity. Hence 
the variational problem (76) yields 



AsJ';0(J'^) > / dsf(r)\A'4,\ 



The energy then satisfies: 
W> dsE|{,|"[K" + J4, 



dr V E ; A^T, 



(80) 



(81) 



Using the definition of k (7f ) and the equilibrium condition 
(65), while writing J,f, in terms of Br and A,f, in terms of Bz 
by inversion of (67), this can be written as 



W > /dsE|{|'{Q| + 



9n 



1 



'■Br dr 



(r'Br) + 



tBz 



/;r5.drj 



(82) 



where Qk = (GM/r^y^^ is the kepler rotation rate and the 
magnetic acceleration g^a is given by (f4). It follows that 
stability is guaranteed if everywhere in the disc we have 



9n 



r±ln{r'Br) 
dr 



r'Bz 



J^rBAr 



(83) 



This may be compared with the estimate (63) from the 
shearing-sheet analysis, which predicts stability against ax- 
isymmetric perturbations if 



9n 



[r^Hr'Br)-\^k,r\^ <Q 



2 



(84) 



The term in kx corresponds to the stabilizing magnetic en- 
ergy term (second term) in (83). While the shearing sheet 
result suggests that the stabilizing effect of the magnetic en- 
ergy term can be removed by taking a sufficiently low wave 
number, the result (83) shows that in fact this term has a fi- 
nite minimum value which can not be removed for any wave 
number. 

As an example consider a disc in which the variables in 
the equilibrium state vary as powers of r, in such a way that 



Bz{r,z = 0) 



(85) 



The potential field corresponding to this Bz has a radial 
component Br{r, z = 0) ~ r"" (cf. Sakurai, f987). Thus, 
the angle of inclination of the field lines at the disc surface 
is constant. For such a field, sufficient condition (83) reduces 
to 



> 0, 



(86) 



independent of p. In equilibrium, the magnetic support gra 
can not exceed the gravity ^k'') that all equilibria which 
vary as powers of r are stable to axisymmetric displace- 
ments. 

To prove stability to axisymmetric perturbations for all 
equilibria in which has constant sign, a stronger version 
of (83) would be needed. This is because the sufficient con- 
dition can be defeated near places where Br has steep gradi- 
ents, by choosing a that peaks near this gradient. As the 
axisymmetric estimate (84) shows, such smaller scale per- 
turbations are stabilized more strongly than suggested by 
(83). 



5 CONCLUSIONS AND DISCUSSION 

A thin disc threaded by a poloidal field which provides part 
of the support against gravity has two kinds of magnetic 
wave modes for motions in the plane of the disc. A com- 
pressive wave is present in which the restoring force is due 
to the potential field outside the disc. It has been studied 
before by ST and Tagger et al. (f990). A radial gradient 
in the field strength introduces modes analogous to internal 
gravity waves. These can become unstable in a sufficiently 
steep gradient in the field strength, and then behave as ex- 
pected from convective or magnetic interchange modes. We 
find that these interchange instabilities in a shearing disc 
behave essentially as expected from an incompressible ap- 
proximation (analogous to the Boussinesq approximation in 
convection) which becomes exact for short wavelengths. 

At longer wavelengths, indications are found for a sec- 
ond form of instability. Its growth maximizes for axisym- 
metric perturbations. In all cases however, it appears only 
when the wavelengths are longer than allowed by the local 
(shearing sheet) approximation under which the results were 
obtained. This form of instability, discussed in some detail in 
Lubow and Spruit (f995), is therefore possibly an artefact. 
A simple global stability bound derived in section 4 also sug- 
gests stability for axisymmetric perturbations. Better global 
stability calculations are needed to settle this issue however. 

In previous work on uniformly rotating discs (ST), the 
instability was found to grow on the rather short magnetic 
time scale gmlr where g^ is the radial acceleration due to 
the magnetic field. The condition for instability was found 
to be that the surface density E and magnetic field strength 
Bz vary such that 



d , E , 



(87) 



The stabilizing effect of the shear flow in the discs stud- 
ied here turns out to be dramatic. We find that significant 
growth of perturbations is possibe only if the growth rate 
in a corresponding system without shear is of the order of 
the dynamical time. The incompressible model gives, as the 
condition for significant growth of linear perturbations: 



Spruit et al: Interchange instability in discs 9 



d , E 



> 25^ 



(88) 



where 5 is the shear rate rdQ/dr. The instability happens 
on short length scales and has algebraic (power law) growth. 
For a disc with finite inner and outer radii, this condition 
also says that the intrinsic (nonrotating) instability time 
scale must be shorter than the shear time scale across the 
whole disc. It is probable that additional globally unstable 
modes would exist in such discs with boundaries. 

Integrations without making the incompressible approx- 
imation agree closely with condition (88). If we assume 
5 Si Q, denote the length scale on which 'E/Bz varies by 
L, and the disc thickness by H, (88) can be written as 



va/cI ^ L/H, 



/? ^ H/L 



(89) 



where va is the Alfven speed in the disc, ~ fi-ff the sound 
speed, and /? the plasma-beta. If 'E/Bz varies smoothly, i.e. 
L ~ r, instability requires the magnetic energy density in- 
side the disc to exceed the geometric mean of the thermal 
and kinetic (orbital) energy densities. If the energy of the ex- 
ternal magnetic field is included, the condition is equivalent 
to the requirement that the magnetic energy of the system 
(disc plus field) must be comparable to the orbital kinetic 
energy. This happens when the magnetic curvature force also 
provides a significant fraction of the support against gravity. 

If, on the other hand, the length scale of the unstable 
gradient is H (the smallest meaningful length scale for this 
gradient), instability requires va<^CsOvP^1. This is 
just the condition at which instability of poloidal fields by 
the Velikhov-Chandrasekhar mechanism (Balbus and Haw- 
ley, 1991) stops operating. Thus it appears that there is a 
significant region in field strengths, 1 ^ va/c., ^ f / H , where 
smoothly varying fields are stable both to interchanges and 
Velikhov-Chandrasekhar-Balbus-Hawley type instabilities. 

Discs in which a poloidal magnetic field provides part 
of the support against gravity are a natural consequence of 
star formation in a magnetic environment. An initially weak 
poloidal field is concentrated by the accretion flow, and its 
energy density will eventually become a significant force op- 
posing accretion, if the field is frozen into the flow. In this 
case, the degree of support is limited by the tendency of 
magnetic fields towards interchange instability, which effec- 
tively allows the field to diffuse radially through the disk. 
Though ambipolar diffusion also has this effect, and may be 
able to leak fields out of the disc during most of its history, 
interchange of some form may well be the main mechanism 
limiting the advection of a poloidal field in the hot inner 
regions of the disc. Our results indicate that quite strong 
poloidal fields may be expected in these regions. Since such 
fields are the main ingredient in the magnetic acceleration 
mechanism of jets and winds, our result lends support to 
the idea that strong poloidal fields are responsible for these 
flows. 
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